#ifndef QMCPLUSPLUS_DATAVIEW_H
#define QMCPLUSPLUS_DATAVIEW_H
#include <OhmmsData/hdf_datatype.h>
#include <mpi/communicator.h>

namespace APPNAMESPACE
{
  template<int D, typename ST=hsize_t>
    struct DataView
    {
      mpi::communicator& mycomm;
      ST count[D];
      ST offset[D];
      ST gcount[D];
      ST stride[D];
      size_t ndata;

      DataView(mpi::communicator& comm)
        : mycomm(comm), ndata(0)
        {
          for(int i=0; i<D; ++i) offset[i]=0;
          for(int i=0; i<D; ++i) stride[i]=1;
        }

      inline size_t size() const { return ndata;}

      /** initialize the data 
       * @param gdims global dimensions
       * @param parts partition
       */
      void init(int *gdims, int *parts, int dir)
      {
        for(int i=0; i<D;++i) gcount[i]=gdims[i];
        ndata=1;
        for(int i=0; i<D;++i) ndata *= count[i] =gdims[i]/parts[i];
        offset[dir]=count[dir]*mycomm.rank();
      }

      /** initialize the data 
       * @param gdims global dimensions
       * @param xpart partition of the zero-th dim
       */
      void init(int *gdims, int xpart)
      {
        for(int i=0; i<D;++i) gcount[i]=gdims[i];
        ndata=count[0]=gdims[0]/xpart;
        for(int i=1; i<D;++i) ndata *= count[i] =gdims[i];
        offset[0]=count[0]*mycomm.rank();
      }

      ////1d partition 
      //void init(int n, bool iscomplex=false)
      //{
      //  count[0]=n/mycomm.size();
      //  offset[0]=count[0]*mycomm.rank();
      //  ndata=count[0];
      //  for(int i=1; i<D; ++i) {count[i]=n; ndata*=n;}
      //  for(int i=0; i<D; ++i) gcount[i]=n;
      //  if(iscomplex) resize_last_dim();
      //}
      ////2d partition 
      //void init(int n, int px, int py, bool iscomplex=false)
      //{
      //  count[0]=n/px; 
      //  count[1]=n/py;
      //  ndata=count[0]*count[1];
      //  offset[0]=count[0]*(mycomm.rank()/py);
      //  offset[1]=count[1]*(mycomm.rank()%py);
      //  for(int i=2; i<D; ++i) {count[i]=n; ndata*=n;}
      //  for(int i=0; i<D; ++i) gcount[i]=n;
      //  if(iscomplex) resize_last_dim();
      //}

      void print(std::ostream& os)
      {
        os << "rank " << mycomm.rank(); 
        os << " offset " ;
        for(int i=0; i<D; ++i) os << " " <<  offset[i];
        os << " count " ;
        for(int i=0; i<D; ++i) os << " " <<  count[i];
        os << " gcount " ;
        for(int i=0; i<D; ++i) os << " " <<  gcount[i];
        os << std::endl;
      }
    };

  template<typename T>
    inline void write_pencils(size_t nitems, size_t nc, T* first, const std::string& rname,
        mpi::communicator& comm)
    {
      char fname[128];
      int ierr;
      MPI_Offset offset=0;
      MPI_Info info;
      MPI_Status status;
      MPI_Datatype etype=mpi::get_mpi_datatype(*first);
      for(int i=0; i<nc; ++i)
      {
        sprintf(fname,"%s.%d",rname.c_str(),i);
        MPI_File fout;
        ierr=MPI_File_open(comm,fname,MPI_MODE_WRONLY|MPI_MODE_CREATE,MPI_INFO_NULL,&fout);
        if(ierr)
        {
          cerr << " Failed to open a mpi file " << fname << endl;
        }
        ierr=MPI_File_set_view(fout,offset,etype,etype,"native",info);
        ierr=MPI_File_write_at_all(fout,offset,(void*)(first+i*nitems),nitems,etype,&status);
        ierr=MPI_File_close(&fout);
      }
    }

  /** write to a hdf5 file
   * @param fname file name
   * @param view data view
   * @param first starting address
   * @param nc number of data sets that share the view
   * @param buffered explicitly set H5Pset buffer
   * @param collective use collective
   */
  template<int D, typename T>
    inline void write(const std::string& rname, const DataView<D>& view, T* first
        , int nc=1, bool buffered=true, bool collective=true)
    {
      char fname[128];
      sprintf(fname,"%s.h5",rname.c_str());
      const hid_t etype=get_h5_datatype(*first);
      MPI_Info info=MPI_INFO_NULL;
      hid_t h_plist = H5Pcreate(H5P_FILE_ACCESS);
      H5Pset_fapl_mpio(h_plist,view.mycomm,info);
      hid_t h_file = H5Fcreate(fname,H5F_ACC_TRUNC,H5P_DEFAULT,h_plist);

      hid_t xfer_plist = H5Pcreate(H5P_DATASET_XFER);
      //Select either collective or indepedent
      if(collective)
        H5Pset_dxpl_mpio(xfer_plist,H5FD_MPIO_COLLECTIVE);
      else
        H5Pset_dxpl_mpio(xfer_plist,H5FD_MPIO_INDEPENDENT);

      //hsize_t buffsize=view.size()/view.count[0];
      hsize_t buffsize=view.size();
      if(buffered)
        hid_t b=H5Pset_buffer(xfer_plist, buffsize, NULL, NULL);

      //hid_t b=H5Pset_chunk (xfer_plist, 1, &buffsize);
      //hid_t b=H5Pset_buffer(xfer_plist, (hsize_t)buffsize, NULL, NULL);
      //unsigned szip_options_mask = H5_SZIP_NN_OPTION_MASK;
      //unsigned szip_pixels_per_block = 16;
      //hid_t b = H5Pset_szip (xfer_plist, szip_options_mask, szip_pixels_per_block);

      hid_t sid1  = H5Screate_simple(D,view.gcount,NULL);
      hid_t memspace=H5Screate_simple(D,view.count,NULL);

      //hsize_t blocks[D];
      //for(int i=0;i<D; ++i) blocks[i]=1;
      herr_t ret;
      for(int i=0; i<nc; ++i)
      {
        sprintf(fname,"u_%d",i);
        hid_t dset_id=H5Dcreate(h_file,fname,etype,sid1,H5P_DEFAULT);
        hid_t filespace=H5Dget_space(dset_id);
        ret=H5Sselect_hyperslab(filespace,H5S_SELECT_SET,view.offset,NULL,view.count,NULL);
        ret = H5Dwrite(dset_id,etype,memspace,filespace,xfer_plist,first+i*view.size());
        H5Sclose(filespace);
        H5Dclose(dset_id);
      }

      H5Sclose(sid1);
      H5Sclose(memspace);
      H5Pclose(xfer_plist);
      H5Pclose(h_plist);
      H5Fclose(h_file);
    }
}
#endif
